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ABSTRACT 

A simple primitive-equation neue for numerical 
prediction of wind fields and temperature fields at 800 and 
400 mb is developed and programmed. Numerical experiments 
are conducted with this model in an effort to suppress the 
Spurious instability encountered in solution of the 
primitive equations. The model proves to be very sensitive 
to spatial finite-differencing schemes. The schemes 
developed in this study tend to suppress the instability 
resulting from truncation and round-off error. Instability 
resulting from unrealistic boundary conditions for a 
limited forecast region is reduced by not allowing 
advection across the boundaries and appears to be con- 
trolled by the finite-differencing schemes used in this 
study. 
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Ne Introduction. 

It is generally recognized that solution of the 
primitive-equations offers the most general and logical 
approach to numerical prediction. Previous attempts at 
solution of the primitive equations by numerical methods 
have met with limited success. However, interest has 
recently been stimulated in the primitive equations by 
unpublished reports indicating that Arakawa [7] has 
succeeded in developing a stable solution. The purpose 
of this study is to develop and program a simple 3-level 
primitive equation model and perform some numerical 
experiments. The main effort in these experiments is 
devoted to the suppression of spurious instability 
encountered in numerical solution of the primitive 
equations. It is hoped that this instability can be 
suppressed by a finite-differencing scheme that will in 


the mean conserve the sum of potential and kinetic energies. 





Fo The Model 

Phillips (1) describes a coordinate system where the 
vertical coordinate p in the x, y, p, t-system is replaced 
by an independent variable VY, where VJ = pa. The following 


relation then holds, moeredns can be x, y, or t: 


(5) - 9 © 9g. 
C/, @& Wak ow (1) 


The subscript p denotes derivatives taken along a constant 
pressure surface. The variable V ranges monotonically from 
zero at the top of the atmosphere to unity at the ground. 

If the above coordinate system is utilized the 
horizontal equation of motion becomes 


Vv 4(v-WV+ TAY +76 -S30 vIr+ 
at aT 1 oT 


r+ ekxvV =o. ©) 


In this study it is assumed that the top of the 
atmosphere P+ is at 200 mb. The atmosphere is divided 


vertically by five levels as follows: 


VV. =0 p=200 mb 
Vi =1/4 p=400 mb 
J, =1/2 p=600 mb 
G=3/h——si—s—s—s—sSsSSSSsSC“‘CNC*CCCO= 8000 mb 
Vy =1 SFC 








For this model Vi» and V3 will be approximated as follows: 


@ @ 
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Vertical velocities are identically zero at the top and 


bottom of the atmosphere (\Vj=0, Vi, =O), reducing vi and 


V3 to 


If the preceding simplifications are utilized and 
equation (2) is applied at the Vand Vi levels, the 


following equations ensue 


ov (ye DY + SA (%-v ) + Vp, ~ 


mae IAT 
Tob V+ FtfTRxy zo (4) 
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SB od; V+ F+rikxeyzo- (5) 
Tr oT 


It is assumed in (4) and (5) that 


OM = 3¥V3 = OW = VE 
oT aT aT AN 
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Moreover the term Via in (4) may be replaced by 
Th 09,6 (23°) Saosin 
eC har, ? ms 


bince We -O andiew —, 
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Wed = OS -F , 
oT AT 


and equation (4) becomes 


Av 49+ Seglts-w) 94 (4 - 


Beir Vit F + KxfEW=O. (6) 
Samy a | 
Similiarly at the Ys level 


G IP = Va. Pa, re Vy Py ~ OF 
oT Ae 


and equation (5) becomes 


a +(¥- V)V + (Ys) + V&t [&, 7 


fe{ Tay - Biba) Job VIT + iF + KXEY = Bi 7: 


Equations (6) and (7) may also be expressed in component 


form as follows: 
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The equation of continuity expressed in Q coordinates 


becomes 
air, V+(7V) + 2 At T) 
at 


This equation can be combined with the thermodynamic 


ee (12) 


equation 


S04 VY VEO+ V 22 = O 
“Te . 
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to yield 


te + vy (TEY) t BAT Sri co. (14) 


Taking the dot product between ]] \Y and equation (2) 


leads to 
BIE OO AV Se my 3 i 
Hag oe ae 
TY VP- TY: Vil =O- (15) 


Multiplying (12) by D and utilizing the relationship 
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ave — WV YO|. Ord? 


results in 
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TV: Vb = VY (ITV) + OST + © 2). 


Substituting the identities 


pat = (os) - Tae. b. ) air 


and 
a : ea PS) ni orerec, : 90D e 
> Are) ger) TN SS 


into (15) and adding (15) to (12) yields a kinetic energy 


equation of the form 








EvV Gg IW y+ Toy) + 3-(1V" r+Tpe)4 
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- 3 (Ov )em — (wa LTY: VITHTTT)S@ = O. (16) 
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The first law of thermodynamics applied to adiabatic 


conditions is 


at at 


Substituting p= TTY eee and expanding the total derivatives 


into Eulerian form for ‘Nj coordinates leads to 


2 (GT) 4+ V(CoT) + TdT - as + 
C aT 3b 


yw V(TT) +: Tatry) 
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Multiplying this equation by 77 and equation (12) by 
Cpl, then adding the resulting equations gives the 


potential energy equation 


d(T) + 7 (70-7 V) + (IT CpT T)- ral adore ; 
Se Ny TC 


eens) + SE =O. (17) 
QgovT 





Combining equations (16) and (17) gives 


d(T T+ a Ve Ve (77 CIV +E V+ TOV) + 
Jt al 


S(TTT+ OV T+ Tht) + Weaesv=0. ae 
Sign ol 26 sei 
Thompson [ 4 | outlines the following three main aspects 


to be considered in solving the primitive equations: 
formulation of boundary conditions, formulation of initial 
conditions and numerical integration of finite difference 
equations. 

The method of numerical integration utilized in this 
Study will be the method of repeated time extrapolation 
outlined by Thompson and shown to be computationally 
Stable for a simple linear model with a grid interval the 
size of the one used in this study if the time interval 
is about 10 min. or less. This method employs a centered 


difference of the form 


Ome Bt 2). At. 


Initial conditions for solution of the primitive 
equations present a real problem. For meteorologically 
significant large-scale motion there is a near balance 
between the coriolis and pressure forces. This means the 
resultant acceleration is an order of magnitude smaller 


than either of these two forces. Small errors in the 
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measurement of initial wind and pressure fields may lead 
to intense divergence fields which, in turn, give large 
pressure changes, thus creating rapid accelerations. 
These spurious accelerations create fictitious gravity 
waves. To exclude these waves initially in this study, 
a form of non-divergent initial winds is used. However, 
this does not necessarily exclude gravity waves for all 
time. 

The forecast region used in this study covers the 
major portion of the Northern Hemisphere. Boundary 
conditions for a limited region such as this are difficult 
to define since the proper type of boundary conditions 
depends on the mathematical character of the primitive 
equations, which have a shifting dual nature, sometimes 
elliptic and sometimes hyperbolic. In this paper all 
parameters are assumed constant with respect to time on 
the boundaries. It is recognized that this assumption 
may lead to the existence of spurious development near the 
boundaries; however, it is felt that if the spatial finite- 
difference schemes are such that in the mean potential and 
kinetic energies are conserved, any boundary instability 
will tend to be averaged out over the forecast region. This 
averaging process will also tend to suppress spurious gravity 
waves resulting from truncation and round-off errors. 

It has been shown by application (see, for example, 


Haltiner [6] ) of Gauss! equation that when an equation 


=~ Ge 





of the same type as (18) is integrated over a closed 
volume the local time change of the sum of kinetic and 
potential energy vanishes. It follows that any finite- 
difference scheme should in the mean conserve the sum of 
the potential and kinetic energies. Specifically, this 
consistency should be achieved separately with respect to 
both vertical and horizontal differencing. 

The vertical consistency is developed as follows. 
First equation (15) is applied at the VY, level. The 


vertical differencing is then performed, leading to 


WEY) ~ V(EW Vth) + LEV vet 


TT, )N T ONE a + |(>,~ O.5)(R + y-977) + 
a (19) 
+b p)TvT | O | 


since YW =0, Vo =0 and V5#V,V,. The derivative ag 
=) 
is approximated by a one-sided difference. A similar 


application of equation (17) results in 


é C ‘ — 
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TG 87 (GEA) ~ SESE - 
rath 


£ (7 +8) B) =O. (20) 
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Adding equations (19) and (20) and comparing the result 
with equation (18) leads to the conclusion that if the 
consistency discussed in the previous paragraph is to be 
maintained, certain terms of equation (19) must cancel 
corresponding terms of equation (20). This inference 
leads to the following equivalences: 


TW + 


x (7+ El By, 


A similar procedure at the V3 level utilizing the 


identity 

v (22) = 3) | To) wee ch 
GT Seika 

leads to 


Vs  & Ty = Dg + O3- Ady 
Tete 


and 


—_ 


k 

= = i = 
>, PD, Cy 207543) T,? r). 
From the previous equations the following relationships can 
be obtained. 

| p.\ k 

T,= (441K) T + 2 (2) T,. 
2 TY, + Fe 1 2 \Veky 3 
Substituting ecco mb, p,=400 mb , Pamoee mb 


k=0.286j/g°k, Nl =% and C,=1.003 j/g°k in the previous system 


=a 





of equations leads to 


2= .643 Ty + .410 T3 (21) 
O, = Pat .442 RTy + .530 RT3 (22) 
Dy = -.058 RT] + .530 RT3 (23) 
Os = Dy+ .058 RT, + .22 RT3 (24) 


If these relationships between geopotential and temperature 
are maintained then the vertical differencing scheme will 
not violate the energy conservation equations. 

With respect to simple derivatives of the form ad 


pee ree in (15), the centered difference approximation 


ao). = Pond De: of 


will obviously vanish when summed over interior points of 
the grid and thus would in the mean conserve potential and 
kinetic energy. However this does not take care of the 
boundaries, nor terms which are not of the simple form 
shown above. 

Since numerical integration is to be performed by 
repeated time extrapolation, the velocity tendencies (8) - 
(11) will have to be computed from oF T and Tx after they 
have been extrapolated forward in time. This will require 
tendency equations for 7/7 and T since D values are expressed 
in terms of T through equations (22) = (24). 


Beecengtion of (12) leads to the approximate form 


ae 2-10-Lriy+\)]. _ 
dC a 


el ge 





From (12) 


Te Seen, 


substituting mit from (25) gives 
o 


w= 40 [ve (ry-ryd . 
The relationship 
“rd = ad 

oT 


utilized in equation (17) yields for the \N, level 
Tepe + Cer | all +9 (ry)/+ Cer V-VT + 
a 


= 


0,-0, ea vive |e sr (P-¢ tGTg)= 0. V0) 
Aims AT 


The bracketed terms of equation (26) may be replaced by the 


following relationship 


st+ Vi (ry) = -7N. 

va aT 

Dividing (26) by IT Cy gives 

207 T G. (T) SO (27) 
2G AY 


A similar approach at the 973 level will result in the 
following tendency equation 

073 + Ve V3 + VG [.68T3 - To - .033T,| =0 (28) 
Jt 3 We : 


A\ 
mae 





based on equations (21) - (24) with assumed zero. 
These tendency equations complete the system required 


for repeated time extrapolation. 
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3. Numerical experiments and results. 
ie Energy Considerations 

In the basic model, the energy balance with respect 
to vertical differencing is maintained through the 
relationship between temperature and geopotential developed 
in Section 2. With respect to horizontal differencing, 
the problem may be considered in two parts, 

A. internal energy balance, 

B. boundary effects. 

For the former problem, partial derivatives of a 
'erosed” form, such as of, can be maintained in balance 
through use of a ey cas differencing scheme. As 
the computation proceeds across the grid, each term is 
added and subtracted once at each grid point, with no 
change in total energy. However, a few derivatives with 
variable coefficients must be evaluated in the horizontal 
acceleration equations, (8) through (11), and also in the 
temperature-change and stability-change equations. In 
these terms, of the form v_J{&, the coefficient varies as 
the grid is traversed, so on the weighting of the term 
at a specific grid point changes between addition and sub- 
traction in central differencing. This can lead to a change 
in total energy across the field, with resultant build- 
up of instability. 


The following approaches were tested on these "open" 
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forms: 

ao lerms left inithas torn. 

This lead to a comparatively short set of forms. For 
example, the bracketed portion of equation (8) was 
simplified to Y19%1 - Wye, - vf. However, this form did 
nocgeswulr 4 1 ee for energy balance, and the 
model became unstable in a very short time with overflow 
occuring in less than four hours. 

b. Terms placed in a "closed" form as a sum of squares. 

The bracketed form of equation (8) now become 

pe 
Vy" -2O(UL+Y1) 
3% 


This form was a distinct improvement on the "open" form; 
however, variable coefficients involved in”, were not 
affected. In any case, there was some instability and this 
model overflowed in about six hours. 

ec. Open form with multipliers replaced by average 
values. 

In this model, the equations of the first approach 
were used, but a "central average" replaced the variable 
coefficient in each case. For example, vdu became, in 


OY 
finite difference notation: 


(vers) tyre) ) Cay tre) 
2 LAY 
This model gave the best results of the three, and was 
used for the remaining experiments. 
The problem of boundary effects proved much more 


difficult. Due to computer storage limitations, this model 
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was programmed using an octagonal grid whose southern 
boundaries are as much as 18° from the equator. Advection 
at these boundaries leads to major problems. 

The following boundary forms were tested: 

a. All time derivatives set equal to zero along the 
boundaries. Since the boundary values were held constant, 
while values at adjacent points were allowed to vary, 
large gradients resulted, giving large-amplitude gravity 
waves which moved inward from all sides at high speed. 
(See Fig. 1) 

b. Boundary values averaged at short time intervals. 

Several smoothing techniques were tried which 
involved only the outside three grid points. Since a 
gravity wave cannot cover the distance between two grid 
points in less than one hour of forecast time, this 
averaging was carried out once each hour. The most 
successful technique involved placing an average of the 
three outer points at the next-to-edge point, then averaging 
this value with that at the outside point and placing the 
result at the outside point. For example; 

(Boundary ) A B C 
becomes SSE aaa? | (A+B+C) /3 C 


This method restricted the instability due to boundary 
effects to the outer two grid points; however, in spite of 
the heavy averaging, large differences developed between 


the outer three grid points after about five hours fore- 


=a 72 





cast time. 

c. All horizontal velocities set equal to zero along 
the boundaries. 

Since the instability developed along the boundaries 
can be traced to energy advected into or out of the forecast 
area, zeroing the horizontal winds at the boundaries, 
thereby zeroing the advection, should solve the problem. 

In practice, it was found that zeroing the horizontal wind 
components at both forecast levels at the outer two grid 
points was effective in controlling boundary-produced araties 
waves. Surprisingly, this did not lead to instability in 

the wind=-component fields themselves in this model. This 
method of boundary condition control was used throughout 


the remaining experiments. 
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Cot rPCcuLon Lemme 

Two friction terms were developed for use with the 
horizontal acceleration equations (8) through (11). For 
horizontal friction, a simple Laplacian multiplied by 10°cmfsec 
was subtracted at each time step. For vertical friction, 
the terms developed by Smagorinsky [| 5] were used, with an 
average skin friction coefficient of 2.2 x 1073 after 
Cressman [ 3]. According to Smagorinsky, the vertical 


friction at any level may be computed from 


vF = ~gdT = g_ (Seti ~ Ses). 
Sp 4p 


Assuming J, =0, 


vF_=g34 ; vWF =g_ (Jy- Ta) 
1h TES 3 AP 


Vertical stress, 7 , is taken as a function of the vertical 


wind shear, so that 


Ret) 


Or Since K is proportional to the magnitude of the wind, 


Mn and € is a constant for any roughness length 


ea ny 


C2 
Q. 


depends on roughness length, and from Cressman's paper 


fm % 


n average value of 2.2 x 1077 was used. Level 4 was taken 
to be the top of the friction level for this term, with 

ivi, estimated as a fraction of the vectorially extrapolated 
wind by the equation 


i 12 Z 2 
in = 0.36 ( (1.192u-0.192u, ) + (1.192v3-0.192v, )”) 


=ao= 





For \y, the cross isobar angle at level 4 was taken to be 
a function of latitude alone by the equation 
ctnd = 1-10° {2f_ y iin sec, 

This gives a cross-isobar angle ranging from near zero 
at the pole to approximately 35° at the base of the grid 
(18N), which is in close agreement with average values. 
The magnitude of , was obtained from the geostrophic 
approximation and the Du, field, which in turn was estimated 
by a simple extrapolation from QD, and o, » namely 

Py =1.192 p, -0.192, , 
Koha So the equation, 

B= ~(PK)2 gh (YW) 


can be readily obtained from the basic J equation. Here 
gp,/OP=1/4z, where z2=5.236 km, the standard depth of the 
800- to 400-mb layer. For (PK),5, the Rossby and 
Montgomery estimate of 500 gm/m sec was used. 

This series of equations was programmed and tested 
on the model. However, the terms involved were so small 
compared to the rest of the model, and the computations 
so time-consuming (approximately 50% increase in computing 


time), that time was not available for a complete evaluation. 


0 
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3. Prognosis of Surface Pressure 

The prognosis of surface pressure was one of the most 
difficult problems encountered because the surface pressure 
change of the basic model was apparently very sensitive to 
small errors in velocity aloft. In fact jy the surface 
pressure field served as a diagnostic chart, Since any 
instability was apparent there long before it appeared 
aloft. 

Four different methods were tried for obtaining the 
prognostic surface pressure field: 

a. Expansion of the basic equation (25). 

The resulting pressure-change field was extremely 
sensitive to changes in velocity divergence at levels l and 
3, and large pressure changes occurred in two to three hours. 

b. Averaging around the grid point. 

In this approach, the products of the pressure value 
and the sum of the wind components at levels one and three 
were differenced across the point being computed. In 


finite difference forn, | 


2 = -4 [Te +) ( U (i+) a Uy (irs) ) ~ TH.) (eaten t Usci-1) } ¥ 


TT a) (Vigne) tigen) - Tow Mg) + VeG-y J | 
This model was not quite as sensitive to divergence 
errors, but still gave unreasonable surface pressure changes 
over periods as short as 6 hours. 


c. Computation of surface pressure directly from the ob 


fielas . 
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<—> } — 





In this case, the deviation of the surface pressure 
from standard was estimated at each time step from the 
deviation of the D fields from standard at 800 and 400 mb 
and the hydrostatic equation; i.e., if b..50 PSoecerimed as 


Viroo - Piers) a etc, then 


Piove = Paoo ye or = OD a os Pave -~ Dgoo 


looo — 
: , 
equation, Ap=Pre, dioos and p=Ppstq tAp. 


/ 
and if we let then from the hydrostatic 
p Prec? zE 
This model was programmed and computed. However, 
Since the p fields are based on temperature and the basic 
inputs to this model were the temperature patterns at 400 
and 800 mb, the resultant surface pressure field lacked 


some of the usual detail. 
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ie Conclusions and suggestions for further study 

In this series of experiments, three finite difference 
forms, three boundary forms, and three surface pressure 
forms were tested. It was found that this model is 
extremely sensitive to finite difference forms used and 
boundary conditions assumed. However, when a new form gave 
a closer approximation to energy balance, stability was 
improved in every case, with the improvement noted being 
in direct proportion to the approach to energy balance. 

It is the opinion of the writers that a stable 
primitive-equation model of this type can be achieved, and 
will eventually become the standard form of numerical 
prognosis. However, faster computers with much larger 
storage capacity must be used to realize this goal. At 
least two more levels of temperature input are needed to 
delineate the geopotential fields ciearly, and the grid 
should be expanded to cover at least a whole hemisphere 
to reduce boundary problems. Even with high speed magnetic 
drum storage, such a model would have a prohibitively long 
running time on a computer such as the CDC 1604; but it 
might be completely feasible on the next generation of 
computers. 

For the present, further experiments with surface 
pressure forms are suggested, possibly using some simpli- 
fied version of the equations currentiy under development 


by Arakawa [ 7] . Also, the friction terms could be 


=2 3= 





streamlined and tested, and a surface terrain term included. 
However, any major revision of this model will require the 
use of a larger, faster computer than any currently 


available. 
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